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In calculating Green functions for interacting quantum systems numerically one often has to resort to finite 
systems which introduces a finite size level spacing. In order to describe the limit of system size going to 
infinity correctly, one has to introduce an artificial broadening larger than the finite size level discretization. 
In this work we compare various discretization schemes for impurity problems, i.e. a small system coupled to 
leads. Starting from a naive linear discretization we will then discuss the logarithmic discretization of the Wilson 
NRG, compare it to damped boundary conditions and arbitrary discretization in energy space. We then discuss 
the importance of choosing the right single particle basis when calculating bulk spectral functions. Finally we 
show the influence of damped boundary conditions on the time evolution of wave packets leading to a NRG- 
tsunami. 



I. INTRODUCTION 

Correlation, or Green, functions are a fundamental concept of condensed matter theory, for an introduction see e.g. Q-dl]. 
However, for interacting systems exact solutions are rare and one often has to resort to perturbative or numerical approaches. In 
calculating Green functions for infinite systems, let it be a bulk Green function, or an impurity problem, numerically one hast 
to resort to a discretized, finite system. In this work we discuss the influence of various discretization schemes on the spectral 
function, i.e. the imaginary part of the retarded Green function. 

To this end we restrict ourselves to non-interacting Fermi systems where numerics can be performed without any approxima- 
tion and the resulting errors can be traced back to the discretization scheme. 

We start with the definition of Green functions in time domain, and derive resolvent equations in frequency domain, which 
demonstrates that for the calculation of spectral functions one does not need the spectrum explicitly. We then discuss the case of 
an impurity problem, namely the resonant level model, where a single level is coupled to a non interacting lead / bath. It turns 
out that it is non-trivial to treat this very simple model accurately on a finite lattice. We then proceed with the problem of an 
energy resolved spectral function of a non-interacting tight binding chain. It turns out, that one can reconstruct the <5-function of 
the spectral function in the continuum. However, care has to be taken in order to avoid discretization errors. Finally we show 
that damped, or Numerical Renormailzation Group- (NRG) like, boundary conditions lead to a phenomenon we call the NRG 
tsunami. 



A. Green Functions in Time Domain 



The lesser (greater) Green functions G < (G > ) and the retarded (advanced) Green functions G r (G a ) are defined d-Hl by 

G>^(M')= -i{A(t)B(t')) (1) 

(2) 



G^(i,t')= -K(B{t')A{t)) 



gv a (m')= -ie(t-o< Mt),B(t') ) = e(t-t')(G> A6 (t,t')-G<^t,t')) 



A.B X 



g a e (t, t') = i e(f - 1)( A(t), B{t') ) = e(f - i)(G< fi (t, f) - g> ~(t, f )) 



(3) 
(4) 



where A and B denote two arbitrary operators. For fermionic operators, ( = —1, [A, B] , = AB + BA denotes the an- 
ticommutator of operators A and B, and for bosonic operators, £ = 1, the anticommutator is replaced by a commutator 
[A, B]_ = AB — BA. Throughout this work (■ • •) denotes the zero temperature ground state average. 



B. Resolvent Representation in Frequency Domain 

In order to simplify notation when switching to frequency domain we assume translational invariance in time and introduce 
the following Green functions of two arbitrary operators A, B: 

G% 6 (t)= -i@(t){A(t)B(0)) (5) 

G^(t)= i@(t){A(0)B(t)), (6) 
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leading to 



G^(t,0)= CGj^-t) - CG^(t) (8) 
G r i,B(*,0)= G+ A (t) + CGT A (t) (9) 

g!b(^o)= CG±,(-t) + G-; .(-*). do) 



The Fourier transformed Green function Q'- I s defined by 



A,B y 

/+oo r+oo 
die-*G+ 5 (t) = -i / dte-*(i(t)S(0)) (11) 
-oo JO 
r+oo 

= -i / dte iaJ * (* |e i:Kt ie- i:Ht i3|*o) (12) 
Jo 

/•+00 

= -i / dt(* |ie i(B °- :K+ " +iT ' ) *5|*o) (13) 
Jo 

B|*o) (14) 



^ Eo-x + v + y, 6 ™ (15) 



and similarly 



(*o|i^^— 5|*o>, (16) 
where a convergence generating 77 = + has been introduced to ensure convergence. Finally, we obtain 

6i,BM = ^M + %H (19) 

^,BM = CGt A (-^) + G^(- w ) (20) 

C. Correction Vector Approach 

It is interesting to note that the resolvent representation Eqs. (TT31 \16[ allows for calculating the Green function without a 
complete knowledge of the spectrum via the correction vector approach |4j]. Starting from a general resolvent expression 

G(E) = (*\A H _ 1 E+ B\V) = (21) 

^ v / 

10 

we obtain the correction vector |£) from the linear system 

(H -E + i V )\0 =m), (22) 

which can be solved by standard solvers. However, note that in this approach one needs a separate run for each desired energy 
E. 
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D. Single Particle Propagator 

The single particle Green functions for fermionic systems are defined by 

G>(x,t;y,t')= -i(c x (t)c+(t')) (23) 

G<(x,t;y,t') = i{c+(t')c x (t)) (24) 

G r (x, t; y, t') = -i Q(t - t'){[c x (t), c+(t')] + ) (25) 

G a (x,t;y,t') = ie(t'-t){[c x (t),4(t')] + ) (26) 

where we use x and y to denote the position in the lattice, c x (t) and c + (t) are the fermionic annihilation and creation operators 
at site x and time t. For spinful calculations x denotes a super index of the spatial coordinate and the spin orbital. From Eqs. ( TT5l 
[Tot we obtain for the retarded and advanced Green functions 

g r (x, y ,u) = g+ . ( w ) - g~+ ( w ) hd 

0°(a;,j/,w) = G r (x,y,Lu)*. (28) 
E. Free Fermions 

Up to now the description for Green functions was completely general. In the following we restrict ourselves to the description 
of non-interacting Fermi systems since our goal is to describe the problem induced by evaluating Green functions on finite 
systems. In result we are able to perform approximation free numerics. Nevertheless, our findings are applicable for interacting 
system and can be exploited by other methods. 

For non-interacting fermions we start with a general Hamiltonian 

Jf = i+-H-d = J2c+H Xty c y . (29) 
y 

We can now switch to a diagonal basis by diagonalizing the matrix H: 

diag(e) = U ■ H ■ U + (30) 
1 = U-U + (31) 
c=U-c, (32) 

where ee are the single particle levels. If the ground state l^o) is non degenerate it is given by 

i*o> = n ^ + i->> < 33 > 

£f<EF 

where ep is the Fermi energy and |— ) is the vacuum state. However, on finite systems at zero temperature this definition of the 
Fermi energy is ambiguous since ep can sit anywhere between the highest occupied and the lowest unoccupied level. We set ef 
in the middle of those two levels to ensure numerical stability. For degenerate ground states one has to take care of the different 
possibilities of filling the highest level. When evaluating expectation values one then has to average degenerate levels at the 
Fermi energy by taking the zero temperature limit of the finite temperature result, e.g. the number of particles N is then given 
by 

N = Urn f GSfo " £f)) Ul x U £iX , (34) 

p— ►oo * — * 

I 

where f() is the fermi function and one should work at a small, but non-vanishing temperature. 

Evaluating the retarded single particle Green functions G r (x, t; y, t 1 ) of Eq. (l27T i in frequency domain using the formulae of 
section UbI we obtain 

g+ H = (V \c x 1 c+|v]/ ) = (1 - (*oN* » UlxUt ' v . (35) 

c *' c v k/Q — K + uj + irj u u — ee + irj 
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and 



g ?c H= -E^ol^l^-^Zl-- (36) 
c v > c x ^— ' W — Si + 1?7 



Finally we obtain from Eqs. d27ll28l ) 



g r {x, y ,u) = g+ + H - g; + ( w ) = £ (37) 

H. RESONANT LEVEL MODEL 

Following the rather general introduction on calculating Green functions we will now concentrate on the spectral function 

A = —%G r {x,x,uj) (38) 

7T 

of a single resonant level edfid coupled via a hybridization t' to a onedimensional lead , where we set the hopping element to 
t = 1. Ignoring the finite width and the cosine dispersion of the lead band results in the wide band limit solution of an area 
normalized Lorentzian 

A=- 9 W , w = t' 2 (39) 

Throughout this section we apply Eq. d37l i to evaluate the resolvent equation d27b . A problem that arises is that the convergence 
factor 77 of d27l i, which is + for continuum leads has to be larger than the finite size level splitting of the leads while it has to be 
much smaller than any physical scale of interest as it also acts as a broadening. In addition, in finite systems the existence of a 
boundary influences the ground state result typically on a scale wbc ~ Ap, where Ap is the level spacing at the Fermi surface. 
Note, that in the case of a single level coupled to a lead with finite width there may exist bound states outside the conduction 
band. While these states are not accessible by single particle of the conduction band due to energy conservation, they can be 
access by few- or many- particle processes. JH. 

The goal of this section is to provide an overview of various discretization schemes of the leads and their impact on the spectral 
function. We start with the natural choice of a finite nearest neighbour tight binding chain. It turns out that their resolution is 
quiet limited. Next we discuss discretization schemes used in the Numerical Renormalization Group (NRG) approach^ and a 
variation of it called Smoothie or DampedJH Boundary conditions, which give a good result for low frequencies, but the high 
frequency results are spoiled. Finally we provide a discretization scheme which is able to provide high resolution on all energy 
scales, where we generalize the variable discretization approach of Nishimoto and Jeckelmann.Jgl 

A. Linear Leads 

In Figure|2]we show the numerical evaluation of the spectral function of a resonant level coupled to a lead via a hybridization 
of t' = 0.1 evaluated for a total number of 700 sites and 350 fermions and hard wall boundary conditions (HWBC). The 

t' t _ t _ t _ * _ t 



Figure 1 : Sketch of a single impurity coupled to a lead via t' and a lead hopping of t. 



corresponding Lorentzian d39l has a half halfwidth of w = t' 2 =0.01 which is already larger than the finite size level splitting 
of Ap sa 0.00352. Nevertheless the figure demonstrates that the discretization of the leads is too coarse to reproduce the 
Lorentzian although the level spacing is already significantly smaller than the resonance width. For an rj larger than the finite 
size level spacing the resonance is artificially broadened while for smaller rj the discrete nature of leads appears in the correlation 
function as can be seen by the spikes of the 77 = 0.0025 result. This observation is explained by looking more closely at the 
resolvent equation d2~7b . In the thermodynamic limit of continuum leads r\ is + and the imaginary part of the resolvent g + + (u>) 

^ Eo-X + u + O+ ^-V-W (40) 
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Figure 2: Spectral function of a resonant level, t' = 0.1, M = 700, N = 350, e d = 0, HWBC, and 77 = 0.05, 0.02, 0.01, 0.0025. L(lu) shows 
the Lorentzian of half halfwidth w = 0.01 corresponding to the spectral function in the wide band limit. 



gives a contribution only at u = E n — Eq, where E n are the eigenenergies of the Hamiltonian "K. By switching to a finite r\ the 
S function gets replaced by a Lorentzian of half halfwidth 77. This corresponds to a convolution of the original spectral function 
with a Lorentzian of the same width and our result is replaced by 

. f°° w/ir 7//-7T 1 w + r] 

LJu)= / de ^ ^ r = 5 ; nr. (41) 

By taking into account the broadening induced by 7] we can fit Eq. (HTT i to our data as shown in Figure[3]and Table[T] The rather 




Figure 3: Lorentzian fits of the spectral function of a resonant level, t' = 0.1, M = 700, N = 350, e d = 0, HWBC, and 77 = 0.05, 0.02, 0.01, 
0.0025. The level spacing at the Fermi energy is Af ~ 0.00352. The full line shows the Lorentzian of half halfwidth w — 0.01 corresponding 
to the spectral function in the wide band limit. 

good results for even large 77 suggest that the results may be strongly improved by unfolding the 77 broadening by a deconvolution. 



B. Logarithmic Discretization 

In order to increase the energy resolution of our lead we can replace the tight binding chain with constant hopping by a tight 
binding chain similar to the one used in NRG where the hopping element is exponentially decreased by a factor A""/ 2 , with 
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0.1 


0.05 


0.02 


0.01 


0.0025 


0.001 


W-L 

l"Log 

WDBC 


0.00942 
0.00855 
0.00945 


0.00973 
0.00876 
0.00974 


0.00990 
0.00887 
0.00987 


0.00995 
0.00891 
0.00979 


0.00999 
0.00894 
0.00942 


0.0170 

0.00895 

0.00929 



Table I: Fits of a Lorentzian J4U to the spectral function of a single impurity coupled to a single lead. All results are in units of the lead 
hopping element t. tol corresponds to the linear leads used in Figure|2]and fits are shown in Figure[3] aiLog corresponds to the logarithmic 
discretization used in Figure|4]and wdbc corresponds to the damped boundary conditions of Figure|6]The analytical result is w = 0.01J. 

A > 1 and n the index of the NRG iteration. |(| Here we use a chain where the hopping is reduced by a factor of A on each bond 
as displayed in Figure[4] 




Figure 4: Sketch of a single impurity coupled to a lead via t' and a NRG like lead hopping of t A n . 

In Figure [5] we show the results for a single level coupled to 3 1 NRG like lead sites where the hopping t = 1 is reduced 
by A -1 = 0.8 on each bond as sketched in Figure |4] For these parameters we get a levelspacing at the Fermi surface of 
Ap ~ 0.000739i. Correspondingly the resolution for tu = is now much higher and the r/ = 0.0025 curve is well resolved 
at small frequencies. However, for large frequencies the resolution drops exponentially leading to spikes even in the r/ = O.Olt 
curve. As a consequence the Lorentzian fits do not work as well as with the (albeit much larger) linear lead, compare Table ©. 
A general feature of the logarithmic discretization consists in the overly broad tails, compare also the the section on frequency 
dependent broadening. 
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Figure 5: Spectral function of a resonant level, = 0, coupled to 31 NRG like lead sites via t' = O.lt using t = 1, A -1 = 0.6. The finite 
size level spacing is Af ~ 0.0007394. L(lj) shows the Lorentzian of half halfwidth w — 0.01 corresponding to the spectral function in the 
wide band limit. 



C. Damped Boundary Conditions 



impurity to Mrs real space sites with constant hopping element t and then coupling to an exponentially decaying lead. In 
FigureQwe show the result for a setup that corresponds to Figure[5]where a hundred additional sites with hopping element t have 
been inserted between the impurity and exponential decaying lead. A similar kind of boundary condition has been originally 
introduced by Vekic and White |7] to mimic the thermodynamic limit in a bulk system. The version employed here was intro- 
duced by Bohr, Schmitteckert, and Wi^ilfle n 1 Ofl to tackle the finite size effects in the evaluation of the Kubo formula for linear 
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Figure 6: Sketch of a single impurity coupled to a lead via t' and a NRG like lead hopping of tA n . 
transport. At a first glance the additional lead sites result in a slight reduction of the spike of the logarithmic discretization only. 
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Figure 7: Spectral function of a resonant level, = 0, coupled to an NRG like Hamiltonian via t' = 0.1 as in Figure[5] A" 1 = 0.8, 
Ma = 30, plus 100 additional sites with a constant hopping of t inserted between the impurity and the NRG like chain leading to a level 
spacing at the Fermi surface of Af « 0.00144. The fits of Lorentzian \4U are shown in Table (U(. The full red line shows the Lorentzian of 
half halfwidth w = 0.01 corresponding to the spectral function in the wide band limit. 

However, one obtains significantly better Lorentzian (HTI) fits compared to the Logarithmic discretization alone. Nevertheless, 
the result is still quite disappointing. 



D. Frequency Dependent Broadening 

In order to remove the spikes in the logarithmic or DBC discretization one has to employ frequency dependent broadening. 
In Figure [8] we present the results for the same system as in Figure [7] However, this time we employ a broadening rj that is 
proportional to the level spacing at energy u. Since the energy spectrum is discrete we take a linear weighted average of the 
level spacing of the level below and above ui. Let us define the level spacing A w at energy us as 

A w = {£„ - e n _i) + (e n+ i - 2e n + e„_i)—— -. r (42) 

0.5 (£„+i - e n _ij 

where n is the level index with 

0.5 (£„_i+£„) < uj < 0.5 (e n + e n+ i) . (43) 

Note, that in the case of degenerate levels one should only count distinct energy levels to avoid a vanishing distance. For an 
w outside the energy range where the corresponding n exists we use the level spacing of the corresponding first or last level 
distance. We then define a relative rj scaling 

t) = * A u . (44) 

For an rj^ = 1 we obtain a broadening 77 that is of the order of the level spacing at energy to. In Figure[8]we compare spectral 
functions for the same system as in Figure|6l only the constant 77 is replaced by a relative 77 broadening. With this approach one 
can eliminate the spikes by using an rj^ > 1.0. However, the tails are still too broad and Eq. (f4TT > can not be used to fit the result, 
as 77 is now energy dependent. In order to filter the 77 induced broadening one would now have to resort to an energy dependent 
deconvolution. 
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Figure 8: Spectral function of a resonant level using DBC as in Figure □ e d = 0, t' = 0.1, A" 1 = 0.8, A/a = 30, Mrs = 30, and 
Ap ~ 0.00144. Here we used a broadening r\ = r\ w * A w , where A u corresponds to the level spacing at energy oj, see Eq. i42\ . The full red 
line shows the Lorentzian of half halfwidth w = 0.01 corresponding to the spectral function in the wide band limit. 



E. Frequency Adapted Grids 



In the previous sections we showed different lattice schemes to evaluate a spectral function. None of the schemes gave actually 
satisfying results. Either the resolution was poor or the tails were not represented correctly. A solution to this problem consists 
in adapting the lattice for each frequency u>. 

We would like to note that this is a generalization to the variable discretization approach of Nishimoto and JeckelmannQ 
in the sense that we can apply a constant broadening for the complete frequency range due to our recipe for constructing 
discretizations. 



1. Momentum Space Leads 



This goal can be achieved by switching to leads in momentum (or energy) space. Let us start with the infinite chain 

M oo 

txCg-Cx-l + ^x^x-l^x ~^ ^x^x^x~l + ^x^x-l^x ' (45) 

x— x— — oo 



We now switch to momentum representations 



oo 

E 

X — — oo 



+ + If* 



(46) 



For a nearest neighbour chain one obtains e(fc) = — 2t cos(fc), however one can now use any desired band e(k). Motivated by 
these considerations we use the general form of a lead in 'momentum space' 

i- [ D \kK{k)e k ~4~c k , (47) 

where N(fc) is the momentum density of states, and D\ ( D u ) the lower (upper) momentum cutoff. Note, here we name 'fc' 
momentum although it can be any labelling. We avoid using an energy density of states since in this work as we might be 
interested in transport properties and keep the flexibility to describe left and right movers by negative and positive momenta. 



Strictly speaking, we should replace the finite lead by a semi-infinite chain. However, this could be incorporated into t x . 
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2. Rediscretization 

In order to use the leads in momentum space for our numerics we have to rediscretize the leads, 

Ml 



Lead 



E 



(48) 



where ki are the Ml discretization points between D\ and D u with fc^_i < kg, which could be used in a scheme as displayed in 
Figure|9] and c e are the fermionic anihilation operators in the new discretization scheme. For convenience we define the interval 
edges 



A I = 

(k e + k e+1 )/2 1 < I < M , 

A, I = M 



(49) 



the level spacing 



(50) 



and density of state weights 



K f .= I dk K(k) . 

Id,-, 



(51) 



In order to preserve the density of states and to ensure canonical commutation relations we have to use the following rule: 
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Figure 9: Sketch of a single impurity coupled to Mrs real space sites which are then coupled to a lead in momentum space representation. 



dk^ik)c(k) -> y — q. 

Z7T Vrf. , V 27T 



(52) 
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which leads to 



2tt f de dk /t—- f d ' r dq 



[c+ c,,] + = -_= / -=V5#) -^qVW)[ct,£ q ] + (53) 



V^" 1 / dk^Wk) dq^qjS^ (54) 



= 5^X7 1 / ' dfcN(fc) (55) 
J dt—i 

= ■ (56) 
Note that in the discretization of a single onedimensional lead we have = 1 and therefore Jit = At. 

3. Level Distribution Function 

In order to obtain a distribution similar to the logarithmic distribution in section HTBl we use an integrated distribution function 
for the levels and discretize the interval [D\, D u ] in equal area sections. Here we use a regularized 1/x function 



p < -W-2 



\J (p-w 2 ) 2 +wj 

Pio S {p, w u w 2 )) = { -j^l ■ -v* <P<™2 (57) 



■ p > w 2 



Plogjp, Wi,W 2 )) 

dkPi os (k,Wi,w 2 )) 



Pi os (p,wi,w 2 ,D h D u )) = B ^ — — . (58) 



to generate the levels. For w 2 = one obtains a logarithmic distribution for the levels e p from P\ og (p, w±, 0, D\, D u )) which 
is slightly smoothed on a scale w\ . The inset of the constant part in the centre enables a linear spacing in the high resolution 
region. This is important for obtaining accurate results in the calculation of linear transport from the Kubo formula. |Q] 

4. RLM Spectral Function 

For a symmetric dispersion e(p) = e(— p) we can now employ another transformation by taking a symmetric level distribution 
for left (p < 0) and right (p > 0) movers and combining them to 

c±, p = (c p ± c_p) 1^2 p>0 (59) 

leading to the lead Hamiltonians 

w± = E e ^) e ±,/±,P' (60) 

where only !K + couples to the impurity. Therefore we can ignore the !K_ part of the Hamiltonian. If we denote the last real 
space site with n, then the coupling to the momentum leads is given by 

-tc+c+ +1 ^-V2tc+Y,V^iC + ,i- (61) 

I 

The main advantage of momentum leads is that one can adapt the discretization to the frequency uj that should be resolved. 
We test this idea by generating a small level distance at frequency u) = —2t cos(p) by slicing P\ og (p, wi, 2wi, D\, D u ) and plot 
the result in Figure [10] Clearly, this approach gives an excellent result which reproduces the central peak and the tails. In 
addition, we can even go back to use a fixed 77, as the level spacing at cj is now always of the same order of magnitude. This 
allows us to fit Eq. ( l4TT i leading to a bare width of wml = 0.00991. In addition, the momentum lead approach allows us to 
restrict the bandwidth to the relevant region. Keeping the parameter of Figure [10] and only changing the momentum cutoff to 
±0.1 we obtain a bare width of wml = 0.00997. 



11 



3 
< 



35 r 
30 - 
25 - 
20 - 
15 - 
10 - 
5 - 





-0.15 



; ~ " ■ " ^-t^a&et 



-0.1 



-0.05 




L(x) 

r|=0.00015 



=ai^^-o o o 



0.05 



0.1 



0.15 



Figure 10: Spectral function of a resonant level using momentum leads and frequency adapted grids. In the numerics we employed a fine 
grained resolution at lu using P\ os (p — w/2, 0.0001, 0.0002, —1, 1), a linearized dispersion eu = 2kt, and a constant broadening r\ = 
0.000154. The momentum leads consist of 50 sites, and the impurity is first coupled to 3 real space sites. The full red line shows the 
Lorentzian of half halfwidth w = 0.01 corresponding to the spectral function in the wide band limit. 



Finally we demonstrate in FigureQT|that within this scheme one can even obtain accurate derivatives of Green functions with 
respect to the frequency. The derivative of ± is given by 



■1 



dui A,B 

d 



(E - Vi + u + iri) 
1 



B|* > 



(62) 
(63) 



It turns out that the numerical evaluation of Eqs. ( 16211631 is more sensitive to the discretization used. We switch to a linear band, 
uj = 2kt, with cutoffs of D\, u = ±0.1, 100 momentum lead sites, the distribution function P\ og (p,wi,2wi), w\ = 0.0005, 
and ?y = 0.0005. As a comparison we plot the exact result and the exact result convoluted with a Lorentzian of half halfwidth 
T) = 0.0005, 



L' v (uj,w) = —L v {uj,w) 



-2uj 



W + 7] 



7T (uj 2 + (w + Tj) 2 ) 



2 ■ 



(64) 



By fitting L' 0005 (u;, w) of Eq. (|4H to the numerical result we obtain w = 0.01004. For the same discretization a fit Eq. (HTI) 
of the spectral function gives w = 0.0099996. 



III. BULK GREEN FUNCTIONS 



In the previous sections we discussed the spectral properties of an impurity, namely a single resonant level, coupled to a non- 
interacting lead. In this case it was natural to look at the local spectral functions of the impurity. We now turn to the momentum 
resolved spectral function of a bulk system. As an example we look at the retarded Green function of a tight binding chain 

% = -f 5^6+gx_i + c+^c, = -2t f dk cos(fc) f+f k 

v X 

At k = ir/2 the spectral function is simply given by A(w) = 5(cu). 
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Figure 11: Derivative of the spectral function of a resonant level using momentum leads. The red line is the analytical result in the wide 
band limit, the blue line is the analytical result convoluted with a Lorentzian of half halfwidth r\ — 0.0005. The circles are obtained from 
momentum leads using a distribution centred around ±cj with Pi og (p, wi, Iwz), w\ = 0.0005, W2 = 0.0001, a linear lead cj(fc) = 2kt with 
cutoffs D\ yU = ±0.1t, and 3 real space sites sandwiched between the impurity and the momentum leads. 



A. Choosing Correct Single Particle States 

For simplicity let us start with a sine function ansatz for the single particle states on M sites 

1 M 

fn = —j= X] sin(A:„a;) c x (65) 

* X — 1 

*» = ^ (66) 

As can be seen in Figure Q~2] this ansatz leads to a fictitious double peak structure and an oscillatory part which does not 
resemble the (5-peak of the system in the thermodynamic limit. The reason for this is that by using the single particle states of 
Eqs. d65ll66*l l we used single particle states which are not eigenstates of the system with periodic boundary conditions (PBC). 

Therefore we repeat this calculation for a tight binding chain employing hard wall boundary conditions (HWBC). The graph 
in Figure [13] shows that the result looks much better now, but we still obtain the fictitious double peak structure. Again, the 
reason for this lies in a wrong single particle basis. However, this time the reason for its failure it is more subtle. The problem is 
that Eq. ( f66b gives the wrong eigenstates for HWBC and the correct single particle basis is given by 



fn = J 1 ^2sin(k n x) c x kn = JJ~fi n=1 ,2,---,M (67) 

x — 1 

leading to the result of Figure [14] which finally resembles the <5-peak broadened by 77. In summary we would like to point out 
that one can obtain nice spectral functions for the bulk system from finite system, however care has to be taken to choose the 
correct representation. Otherwise spurious structures may appear. The need for using the sine solution for HWBC has been 
pointed out by Benthien, Gebhard, and Jeckelmann lfTTll . It was shown by Ulbricht and Schmitteckert that for interacting particle 
in a harmonic trap one can obtain spectral functions from finite system by resorting to hermite polynomials. 1 1211 



IV. POOR MAN'S DECONVOLUTION: SELF ENERGY SHARPENING 



In this section we discuss a simple strategy to "sharpen" Green functions obtained in the previous sections. It was first applied 
to improve the spectral function of a spin polarized ondedimensional Hubbard model. HH As an example we improve the result 
for the spectral function of the tight binding chain of the previous section. A straight forward strategy to improve the result 
consists in a deconvolution to remove the broadening introduced by the finite 7/ in the denominator of the resolvents, i.e. perform 
the inverse operation of Eq. ( RTV While this procedure is mathematically well defined, it is highly unstable numerically and very 
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Figure 12: Spectral function of the tight binding chain at k = tt/2 using the single particle states of Eqs. 
with M — 150 lattice sites and periodic boundary conditions (PBC) and sine function as the basis set. 
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sensitive to numerical noise. 1 13H Here we want to introduce a simple method to improve the results, which has the advantage of 
being stable and allows to preserve analytical properties of spectral functions, e.g. A(u>) < 0. 
We start by defining the self energy for the infinite system 

0» = V A , - n+ (68) 

u> — zj(wJ + i0+ 

and the self energy S r; (w) for the finite system which includes the finite broadening rj 

1 



w - T, v {w) + iCH 



(69) 



When switching from the infinite system to the finite system we should trace out the discarded part which would lead to a 
contribution to the self energy of the reduced system. However, we do not know this part. Therefore we make the ansatz that 
the discarded part can be modeled by the vq term we already have in the resolvents to enable the mixing of excited states. In this 
way we obtain our "poor man's deconvolution" ansatz S r) = £ + vq: 

1 



Y, v {uj) 
E(w) 



1 



i0^ 



irj + i0^ 



(70) 



(71) 



where G~(w) is the Green function we obtain from our numerics on a finite lattice and rj is the broadening explicitly used in the 
numerics. The result of this sharpening applied to the data of Figure[l4]is shown in Figure[l5] Since real and imaginary part 
are vanishing up to numerical precision we have rediscovered the <5-peak structure of the infinite system from the numerics on a 
finite lattice. 



V. THE NRG TSUNAMI 



Finally we would like to discuss the effect of damped boundary conditions on the dynamics of wave packets to show that 
boundary conditions are not only important in frequency space, but that they can also change the results in time domain dramat- 
ically. In this section we follow the route of [ 14j] adapted for free fermions using exact diagnalization of the quadratic form. In 
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Figure 13: Spectral function of the tight binding chain at k = tt/2 using the single particle states of Eqs. J65U66t and a tight binding chain 
with M — 150 lattice sites and periodic boundary conditions (PBC) and sine function as the basis set. 



Figure[T6]we show the initial states of a 201 site tight binding chain with periodic boundary conditions at half filling, where half 
a fermion was trapped in the middle of the system by applying a Gaussian potential with a width of a = 2.5. In addition we 
show the response to the same perturbation for a system where the PBC are replaced by damped boundary conditions, compare 
Figure [6] 

In Figure [17] we show the system after performing a time evolution up to time T = 90, in which the homogeneous system 
moves a distance of T * vp = 180 sites. Due to the periodic boundary conditions the wave packets appear now at position x « 
100 ± 20. Note that by applying a low energy perturbation we create excitations around ±kp and our initial wave packet splits 
into a right and left moving wave packet travelling at ±vp, where kp is the Fermi momentum and vp is the excitation velocity. 
However, once the wave packet hits the region of damped boundary conditions (DBC) it sees an exponentially decreasing 
hopping element which results in an exponentially decreasing excitation velocity. In the similar way as long water waves hitting 
the shore, where the excitation velocity is decreased, our wave packet has to pile up like a tsunami, since the front is moving 
much slower than the back. In addition each changed hopping element creates a small back reflection leading to an additional 
wiggling. 



VI. SUMMARY 



In summary we have shown that one can extract sound results for the thermodynamic limit in the context of calculating 
Green functions from a finite lattice. However, care has to be taken to choose the correct representation of states. In order to 
calculate spectral functions at finite frequencies for impurity problems one should apply frequency adapted grids to achieve high 
resolution. A sharpening procedure based on the self energy can be used as a replacement for deconvoluting the broadening 
introduced by a finite r}. 
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Figure 14: Spectral function of the tight binding chain at k = n/2 using the single particle states of Eq. J67t and a tight binding chain with 
M = 150 lattice sites and periodic boundary conditions (PBC) and sine function as the basis set. 
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Figure 16: Density of a M — 201 site tight binding chain with N = 101 fermions, where a Gaussian potential of width a = 2.5 was applied to 
trap an additional half fermion in the centre region in addition to an average density of p = 0.5. In the results for periodic boundary conditions 
(PBC) are given by the crosses. The result for Damped Boundary conditions as displayed in Figure [6] applied to the left and right end of the 
chain is displayed by the line with plusses. 
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Figure 17: Density of a M — 201 site tight binding chain with TV = 101 fermions, where a Gaussian potential of width a = 2.5 as in 
Figure [T6l after performing a time evolution of time T = 90 with the Hamiltonian without the perturbation. The results for periodic boundary 
conditions (PBC) are given by the crosses. The result for damped boundary conditions as displayed in Figure [6] applied to the left and right 
end of the chain is shown by the line with plusses. 



